1 Setup

1.1 Load libraries

library(survival)
library(piecewiseSEM)
source("spat1f/init_analysis.R")
Data last updated on 2024-05-03

1.2 Handy functions not in spat1f

trend_by_pre_wt <- function(model, term){
  emmeans::emtrends(model, var = "cat_pre_wt_log_scale", specs = term) %>% 
  pairs()
}

contrast_by_pre_wt <- function(model, term, type = c("pairwise", "trt.vs.ctrl"), at = c(-2, 2)){
  type <- match.arg(type)
  f <- as.formula(sprintf("%s ~ %s | cat_pre_wt_log_scale", type, term))
  emmeans::emmeans(model, 
                 f, 
                 var = term, 
                 at = list(cat_pre_wt_log_scale = at))
}

check_mod <- function(model){
 DHARMa::simulateResiduals(model) %>% 
    plot()
  performance::posterior_predictive_check(model)
}

1.3 Quick cleaning

# Main data.frame for analyses
d <- ref_data %>% 
  filter(error == 0) %>% 
  mutate(
    cat_pre_wt_log = log(cat_pre_wt),
    cat_pre_wt_log_scale = as.numeric(scale(log(cat_pre_wt))),
    cat_size = ifelse(cat_pre_wt < median(cat_pre_wt), "small", "big"),
    has_var = ifelse(var_trt != "constant", 1, 0),
    mean_toxic_conc_scale = as.numeric(scale(mean_toxic_conc)),
    mean_toxic_scale = as.numeric(scale(mean_toxic)),
    area_herb_log_scale = as.numeric(scale(log(area_herb+1))),
    var_toxic_12_scale = as.numeric(scale(var_toxic_12)),
    on_toxic_logit_scale = as.numeric(scale(adjust_prop(on_toxic, 
                                                        trans = "emp", 
                                                        nudge.method = "none", 
                                                        na.action = "ignore"))),
    ava_qual_scale = as.numeric(scale(ava_qual)),
    ava_qual_logit_scale = as.numeric(scale(adjust_prop(ava_qual, 
                                                        trans = "emp", 
                                                        nudge.method = "none", 
                                                        na.action = "ignore"))),
    sl_mean_obs_log_scale = as.numeric(scale(log(sl_mean_obs))),
    cat_pre_wt_log_scale_sq = as.numeric(scale(log(cat_pre_wt)^2)),
    RGR_scale = as.numeric(scale(RGR)), 
    sl_mean_pred1 = shape_1 * scale_1,
    sl_mean_pred2 = shape_2 * scale_2,
    sl_mean_pred3 = shape_3 * scale_3,
    sl_mean_pred4 = shape_4 * scale_4,
    k1 = scale_1 / scale_2,
    k2 = scale_3 / scale_4,
    prop_explore_logit = qlogis(prop_explore),
    prop_explore_logit_scale = as.numeric(scale(prop_explore_logit))
  )

# Subset of data.frame for SEM
d2 <- d %>% 
  filter(
    var_trt != "constant"
  ) %>% 
  filter(
    !is.na(area_herb_log_scale) & 
      !is.na(var_toxic_12_scale) & 
      !is.na(mean_toxic_conc_scale) &
      !is.na(sl_mean_obs) &
      !is.na(prop_explore) & 
      !is.na(RGR)) %>% 
  mutate(
    var_high = ifelse(var_trt == "high_var", 1, 0),
    beta_red = ifelse(as.character(beta) == 5, 1, 0),
    beta_white = ifelse(as.character(beta) == 0, 1, 0),
    beta_red_cat_pre_wt_log_scale = beta_red * cat_pre_wt_log_scale,
    beta_white_cat_pre_wt_log_scale = beta_white * cat_pre_wt_log_scale,
    var_high_cat_pre_wt_log_scale = var_high * cat_pre_wt_log_scale,
    beta_numeric_scale = as.numeric(scale(as.numeric(beta)))
  )

# Subset of data.frame for more detailed analyses
d3 <- d %>% 
  filter(
    var_trt != "constant"
  ) %>% 
  filter(
    sl_mean_pred2 < 2000 &
    sl_mean_pred1 > 0 & 
      sl_mean_pred2 > 0 &
      sl_mean_pred3 > 0 & 
      sl_mean_pred4 > 0 &
      scale_1 > 0 & 
      shape_1 > 0 & 
      scale_2 > 0 & 
      shape_2 > 0 & 
      scale_3 > 0 & 
      shape_3 > 0 & 
      scale_4 > 0 & 
      shape_4 > 0 & 
      !is.na(area_herb_log_scale) & 
      !is.na(var_toxic_12_scale) & 
      !is.na(mean_toxic_conc_scale) &
      !is.na(on_toxic_logit_scale) &
      !is.na(ava_qual_scale) &
      !is.na(sl_mean_obs) &
      !is.na(prop_explore) & 
      !is.na(RGR)) %>% 
  mutate(
    var_high = ifelse(var_trt == "high_var", 1, 0),
    beta_red = ifelse(as.character(beta) == 5, 1, 0),
    beta_white = ifelse(as.character(beta) == 0, 1, 0),
    beta_red_cat_pre_wt_log_scale = beta_red * cat_pre_wt_log_scale,
    beta_white_cat_pre_wt_log_scale = beta_white * cat_pre_wt_log_scale,
    var_high_cat_pre_wt_log_scale = var_high * cat_pre_wt_log_scale,
    beta_numeric_scale = as.numeric(scale(as.numeric(beta)))
  )

2 Herbivore Performance

2.1 RGR

2.1.1 Model summary

rgr_m <- glmmTMB(
  RGR ~ 
    (var_trt + beta) * cat_pre_wt_log_scale +
    I(cat_pre_wt_log_scale^2) +  # residual plot show significant non-linearity
    (1|session_id),
  data = d %>% 
    filter(
      var_trt != "constant"
    ) 
); summary(rgr_m)
 Family: gaussian  ( identity )
Formula:          
RGR ~ (var_trt + beta) * cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) +  
    (1 | session_id)
Data: d %>% filter(var_trt != "constant")

     AIC      BIC   logLik deviance df.resid 
  -879.2   -849.4    450.6   -901.2      100 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev.
 session_id (Intercept) 3.979e-06 0.001995
 Residual               1.603e-05 0.004004
Number of obs: 111, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 1.6e-05 

Conditional model:
                                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)                          0.0162504  0.0013457  12.075  < 2e-16 ***
var_trtlow_var                      -0.0034872  0.0007740  -4.506 6.62e-06 ***
beta0                               -0.0003065  0.0009664  -0.317  0.75112    
beta5                               -0.0013545  0.0009173  -1.477  0.13978    
cat_pre_wt_log_scale                -0.0048759  0.0011991  -4.066 4.78e-05 ***
I(cat_pre_wt_log_scale^2)           -0.0021953  0.0006731  -3.262  0.00111 ** 
var_trtlow_var:cat_pre_wt_log_scale  0.0025873  0.0007897   3.276  0.00105 ** 
beta0:cat_pre_wt_log_scale           0.0011924  0.0009491   1.256  0.20897    
beta5:cat_pre_wt_log_scale           0.0024731  0.0009710   2.547  0.01087 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(rgr_m, type = "III")
Analysis of Deviance Table (Type III Wald chisquare tests)

Response: RGR
                                Chisq Df Pr(>Chisq)    
(Intercept)                  145.8141  1  < 2.2e-16 ***
var_trt                       20.3008  1  6.617e-06 ***
beta                           2.3963  2   0.301758    
cat_pre_wt_log_scale          16.5348  1  4.777e-05 ***
I(cat_pre_wt_log_scale^2)     10.6374  1   0.001108 ** 
var_trt:cat_pre_wt_log_scale  10.7334  1   0.001052 ** 
beta:cat_pre_wt_log_scale      6.4936  2   0.038899 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.1.2 Post-hoc contrasts

trend_by_pre_wt(rgr_m, "var_trt")
 contrast           estimate      SE  df t.ratio p.value
 high_var - low_var -0.00259 0.00079 100  -3.276  0.0014

Results are averaged over the levels of: beta 
contrast_by_pre_wt(rgr_m, "var_trt")
$emmeans
cat_pre_wt_log_scale = -2:
 var_trt     emmean      SE  df  lower.CL upper.CL
 high_var  0.014224 0.00355 100  0.007178  0.02127
 low_var   0.005562 0.00312 100 -0.000635  0.01176

cat_pre_wt_log_scale =  2:
 var_trt     emmean      SE  df  lower.CL upper.CL
 high_var -0.000393 0.00208 100 -0.004512  0.00373
 low_var   0.001295 0.00238 100 -0.003421  0.00601

Results are averaged over the levels of: beta 
Confidence level used: 0.95 

$contrasts
cat_pre_wt_log_scale = -2:
 contrast           estimate      SE  df t.ratio p.value
 high_var - low_var  0.00866 0.00181 100   4.790  <.0001

cat_pre_wt_log_scale =  2:
 contrast           estimate      SE  df t.ratio p.value
 high_var - low_var -0.00169 0.00171 100  -0.988  0.3256

Results are averaged over the levels of: beta 
trend_by_pre_wt(rgr_m, "beta")
 contrast         estimate       SE  df t.ratio p.value
 (beta-5) - beta0 -0.00119 0.000949 100  -1.256  0.4232
 (beta-5) - beta5 -0.00247 0.000971 100  -2.547  0.0329
 beta0 - beta5    -0.00128 0.000940 100  -1.363  0.3644

Results are averaged over the levels of: var_trt 
P value adjustment: tukey method for comparing a family of 3 estimates 
contrast_by_pre_wt(rgr_m, "beta")
$emmeans
cat_pre_wt_log_scale = -2:
 beta    emmean      SE  df  lower.CL upper.CL
 -5    0.012890 0.00382 100  0.005313  0.02047
 0     0.010199 0.00333 100  0.003591  0.01681
 5     0.006589 0.00319 100  0.000255  0.01292

cat_pre_wt_log_scale =  2:
 beta    emmean      SE  df  lower.CL upper.CL
 -5   -0.001439 0.00226 100 -0.005922  0.00304
 0     0.000639 0.00231 100 -0.003940  0.00522
 5     0.002153 0.00257 100 -0.002946  0.00725

Results are averaged over the levels of: var_trt 
Confidence level used: 0.95 

$contrasts
cat_pre_wt_log_scale = -2:
 contrast         estimate      SE  df t.ratio p.value
 (beta-5) - beta0  0.00269 0.00221 100   1.218  0.4456
 (beta-5) - beta5  0.00630 0.00223 100   2.820  0.0158
 beta0 - beta5     0.00361 0.00211 100   1.711  0.2062

cat_pre_wt_log_scale =  2:
 contrast         estimate      SE  df t.ratio p.value
 (beta-5) - beta0 -0.00208 0.00205 100  -1.015  0.5689
 (beta-5) - beta5 -0.00359 0.00206 100  -1.746  0.1935
 beta0 - beta5    -0.00151 0.00211 100  -0.719  0.7529

Results are averaged over the levels of: var_trt 
P value adjustment: tukey method for comparing a family of 3 estimates 

2.1.3 Misc stats

r2(rgr_m, tolerance = 10^-10)
# R2 for Mixed Models

  Conditional R2: 0.555
     Marginal R2: 0.444

2.1.4 Check model

check_mod(rgr_m)

2.1.5 Plot

g1 <- marginal_effects(rgr_m, terms = c("cat_pre_wt_log_scale","beta")) %>% 
  ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = yhat)) + 
  geom_ribbon(aes(ymax = upper, ymin = lower, fill = beta), alpha = 0.2) + 
  geom_line(aes(color = beta), linewidth = 2) + 
  geom_point(
    data = rgr_m$frame,
    aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = RGR, color = beta),
    size = 3, alpha = 0.8
  ) + 
  theme_bw(base_size = 15) + 
  theme(legend.position = "top") +
  scale_x_continuous(trans = "log10", labels = fancy_scientific) +
  scale_color_brewer(type = "qual", aesthetics = c("fill","color")) + 
  labs(x = "Cat pre-weight (g)", y = expression(RGR~(hour^-1)),
       color = expression(beta), fill = expression(beta))

g2 <- marginal_effects(rgr_m, terms = c("cat_pre_wt_log_scale","var_trt")) %>% 
  ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = yhat)) + 
  geom_ribbon(aes(ymax = upper, ymin = lower, fill = var_trt), alpha = 0.2) + 
  geom_line(aes(color = var_trt), linewidth = 2) + 
  geom_point(
    data = rgr_m$frame,
    aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = RGR, color = var_trt),
    size = 3, alpha = 0.8
  ) + 
  theme_bw(base_size = 15) + 
  scale_x_continuous(trans = "log10", labels = fancy_scientific) +
  theme(legend.position = "top") +
  scale_color_discrete(type = c("navy","skyblue"), 
                       aesthetics = c("fill","color"), 
                       label = c("High", "Low")) + 
  labs(x = "Cat pre-weight (g)", y = expression(RGR~(hour^-1)),
       color = "Variation", fill = "Variation")


ggarrange(g1, g2 + labs(y = ""))

2.2 Time to pupation

2.2.1 Model summary

pupt_m_full <- glmmTMB(
  time_to_pupation ~ 
    (var_trt + beta) * cat_pre_wt_log_scale + 
    I(cat_pre_wt_log_scale^2) + # Residuals show significant non-linearity
    (1|session_id),
  data = d %>% 
    filter(var_trt != "constant"), 
  family = Gamma(link = "log"),
); summary(pupt_m_full)
 Family: Gamma  ( log )
Formula:          time_to_pupation ~ (var_trt + beta) * cat_pre_wt_log_scale +  
    I(cat_pre_wt_log_scale^2) + (1 | session_id)
Data: d %>% filter(var_trt != "constant")

     AIC      BIC   logLik deviance df.resid 
   338.1    366.3   -158.0    316.1       85 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 2.839e-11 5.328e-06
Number of obs: 96, groups:  session_id, 5

Dispersion estimate for Gamma family (sigma^2): 0.0254 

Conditional model:
                                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)                          2.180302   0.036853   59.16  < 2e-16 ***
var_trtlow_var                       0.102145   0.033998    3.00  0.00266 ** 
beta0                               -0.001845   0.040823   -0.05  0.96395    
beta5                                0.059543   0.041181    1.45  0.14821    
cat_pre_wt_log_scale                -0.428662   0.032643  -13.13  < 2e-16 ***
I(cat_pre_wt_log_scale^2)           -0.081267   0.019069   -4.26 2.03e-05 ***
var_trtlow_var:cat_pre_wt_log_scale  0.026518   0.038876    0.68  0.49516    
beta0:cat_pre_wt_log_scale          -0.030363   0.043195   -0.70  0.48210    
beta5:cat_pre_wt_log_scale          -0.126801   0.047095   -2.69  0.00709 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(pupt_m_full, type = "III")
Analysis of Deviance Table (Type III Wald chisquare tests)

Response: time_to_pupation
                                 Chisq Df Pr(>Chisq)    
(Intercept)                  3500.2207  1  < 2.2e-16 ***
var_trt                         9.0264  1   0.002661 ** 
beta                            2.7351  2   0.254736    
cat_pre_wt_log_scale          172.4448  1  < 2.2e-16 ***
I(cat_pre_wt_log_scale^2)      18.1614  1  2.029e-05 ***
var_trt:cat_pre_wt_log_scale    0.4653  1   0.495162    
beta:cat_pre_wt_log_scale       7.5589  2   0.022836 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pupt_m <- glmmTMB(
  time_to_pupation ~ 
    var_trt + beta * cat_pre_wt_log_scale + 
    I(cat_pre_wt_log_scale^2) + # Residuals show significant non-linearity
    (1|session_id),
  data = d %>% 
    filter(var_trt != "constant"), 
  family = Gamma(link = "log"),
); summary(pupt_m)
 Family: Gamma  ( log )
Formula:          
time_to_pupation ~ var_trt + beta * cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) +  
    (1 | session_id)
Data: d %>% filter(var_trt != "constant")

     AIC      BIC   logLik deviance df.resid 
   336.5    362.2   -158.3    316.5       86 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 2.159e-11 4.647e-06
Number of obs: 96, groups:  session_id, 5

Dispersion estimate for Gamma family (sigma^2): 0.0255 

Conditional model:
                            Estimate Std. Error z value Pr(>|z|)    
(Intercept)                 2.180555   0.036937   59.03  < 2e-16 ***
var_trtlow_var              0.107447   0.033197    3.24  0.00121 ** 
beta0                      -0.003578   0.040842   -0.09  0.93019    
beta5                       0.059536   0.041271    1.44  0.14915    
cat_pre_wt_log_scale       -0.420147   0.030221  -13.90  < 2e-16 ***
I(cat_pre_wt_log_scale^2)  -0.083365   0.018881   -4.42 1.01e-05 ***
beta0:cat_pre_wt_log_scale -0.026011   0.042849   -0.61  0.54383    
beta5:cat_pre_wt_log_scale -0.126834   0.047183   -2.69  0.00719 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(pupt_m, type = "III")
Analysis of Deviance Table (Type III Wald chisquare tests)

Response: time_to_pupation
                              Chisq Df Pr(>Chisq)    
(Intercept)               3485.0791  1  < 2.2e-16 ***
var_trt                     10.4759  1   0.001209 ** 
beta                         2.8053  2   0.245947    
cat_pre_wt_log_scale       193.2743  1  < 2.2e-16 ***
I(cat_pre_wt_log_scale^2)   19.4946  1  1.009e-05 ***
beta:cat_pre_wt_log_scale    7.6890  2   0.021397 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.2.2 Post-hoc contrasts

contrast_by_pre_wt(pupt_m, "var_trt")
$emmeans
cat_pre_wt_log_scale = -2:
 var_trt  emmean     SE  df asymp.LCL asymp.UCL
 high_var  2.808 0.0769 Inf     2.657      2.96
 low_var   2.915 0.0818 Inf     2.755      3.08

cat_pre_wt_log_scale =  2:
 var_trt  emmean     SE  df asymp.LCL asymp.UCL
 high_var  0.924 0.0672 Inf     0.792      1.06
 low_var   1.031 0.0725 Inf     0.889      1.17

Results are averaged over the levels of: beta 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
cat_pre_wt_log_scale = -2:
 contrast           estimate     SE  df z.ratio p.value
 high_var - low_var   -0.107 0.0332 Inf  -3.237  0.0012

cat_pre_wt_log_scale =  2:
 contrast           estimate     SE  df z.ratio p.value
 high_var - low_var   -0.107 0.0332 Inf  -3.237  0.0012

Results are averaged over the levels of: beta 
Results are given on the log (not the response) scale. 
trend_by_pre_wt(pupt_m, "beta")
 contrast         estimate     SE  df z.ratio p.value
 (beta-5) - beta0    0.026 0.0428 Inf   0.607  0.8163
 (beta-5) - beta5    0.127 0.0472 Inf   2.688  0.0197
 beta0 - beta5       0.101 0.0473 Inf   2.132  0.0834

Results are averaged over the levels of: var_trt 
P value adjustment: tukey method for comparing a family of 3 estimates 
contrast_by_pre_wt(pupt_m, "beta")
$emmeans
cat_pre_wt_log_scale = -2:
 beta emmean     SE  df asymp.LCL asymp.UCL
 -5    2.741 0.0941 Inf     2.557      2.93
 0     2.790 0.0971 Inf     2.599      2.98
 5     3.054 0.1088 Inf     2.841      3.27

cat_pre_wt_log_scale =  2:
 beta emmean     SE  df asymp.LCL asymp.UCL
 -5    1.061 0.0845 Inf     0.895      1.23
 0     1.005 0.0805 Inf     0.847      1.16
 5     0.866 0.0917 Inf     0.687      1.05

Results are averaged over the levels of: var_trt 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
cat_pre_wt_log_scale = -2:
 contrast         estimate     SE  df z.ratio p.value
 (beta-5) - beta0  -0.0484 0.1018 Inf  -0.476  0.8828
 (beta-5) - beta5  -0.3132 0.1133 Inf  -2.763  0.0158
 beta0 - beta5     -0.2648 0.1136 Inf  -2.331  0.0516

cat_pre_wt_log_scale =  2:
 contrast         estimate     SE  df z.ratio p.value
 (beta-5) - beta0   0.0556 0.0875 Inf   0.636  0.8005
 (beta-5) - beta5   0.1941 0.0915 Inf   2.122  0.0854
 beta0 - beta5      0.1385 0.0928 Inf   1.492  0.2945

Results are averaged over the levels of: var_trt 
Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 

2.2.3 Misc stats

r2(pupt_m, tolerance = 10^-10)
Random effect variances not available. Returned R2 does not account for random effects.
# R2 for Mixed Models

  Conditional R2: NA
     Marginal R2: 0.881

2.2.4 Check model

check_mod(pupt_m)

2.2.5 Plot

g3 <- marginal_effects(pupt_m, terms = c("cat_pre_wt_log_scale","beta")) %>% 
  ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = yhat)) + 
  geom_ribbon(aes(ymax = upper, ymin = lower, fill = beta), alpha = 0.2) + 
  geom_line(aes(color = beta), linewidth = 2) + 
  geom_point(
    data = pupt_m$frame,
    aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = time_to_pupation, color = beta),
    size = 3, alpha = 0.8
  ) + 
  theme_bw(base_size = 15) + 
  theme(legend.position = "top") +
  scale_y_continuous(trans = "log2") +
  scale_x_continuous(trans = "log10", labels = fancy_scientific) +
  scale_color_brewer(type = "qual", aesthetics = c("fill","color")) + 
  labs(x = "Cat pre-weight (g)", y = "Time to pupation (days)",
       color = expression(beta), fill = expression(beta))

g4 <- marginal_effects(pupt_m, terms = c("var_trt")) %>% 
  ggplot(aes(x = var_trt, y = yhat)) + 
  geom_point(
    data = pupt_m$frame,
    aes(x = var_trt, y = time_to_pupation, color = var_trt),
    position = position_jitter(width = 0.2),
    size = 3, alpha = 0.8
  ) + 
  geom_pointrange(
    aes(ymax = upper, ymin = lower), 
    color = "black",
    size = 1, linewidth = 2, shape = 1,
  ) +
  theme_bw(base_size = 15) + 
  scale_y_continuous(trans = "log2") +
  theme(legend.position = "top") +
  scale_color_discrete(type = c("navy","skyblue"), 
                       aesthetics = c("color"), 
                       label = c("High", "Low")) + 
  labs(x = "Variation", y = "Time to pupation (days)",
       color = "Variation")


ggarrange(g3, g4 + labs(y = ""))

2.3 Eclosure

2.3.1 Model summary

eclosure_m_full <- glmmTMB(
  eclosed ~ 
    (var_trt + beta) * cat_pre_wt_log_scale + 
    (1|session_id),
  data = d %>% 
    filter(var_trt != "constant"), 
  family = binomial(link = "logit"),
); summary(eclosure_m_full)
 Family: binomial  ( logit )
Formula:          
eclosed ~ (var_trt + beta) * cat_pre_wt_log_scale + (1 | session_id)
Data: d %>% filter(var_trt != "constant")

     AIC      BIC   logLik deviance df.resid 
   185.3    210.9    -83.6    167.3      119 

Random effects:

Conditional model:
 Groups     Name        Variance Std.Dev.
 session_id (Intercept) 0.1661   0.4076  
Number of obs: 128, groups:  session_id, 5

Conditional model:
                                    Estimate Std. Error z value Pr(>|z|)
(Intercept)                          0.53530    0.41050   1.304    0.192
var_trtlow_var                       0.01707    0.37507   0.046    0.964
beta0                               -0.21544    0.46265  -0.466    0.641
beta5                               -0.50749    0.45596  -1.113    0.266
cat_pre_wt_log_scale                 0.10452    0.36089   0.290    0.772
var_trtlow_var:cat_pre_wt_log_scale  0.17526    0.36489   0.480    0.631
beta0:cat_pre_wt_log_scale           0.24778    0.43553   0.569    0.569
beta5:cat_pre_wt_log_scale           0.20159    0.44699   0.451    0.652
car::Anova(eclosure_m_full, type = "III")
Analysis of Deviance Table (Type III Wald chisquare tests)

Response: eclosed
                              Chisq Df Pr(>Chisq)
(Intercept)                  1.7004  1     0.1922
var_trt                      0.0021  1     0.9637
beta                         1.2502  2     0.5352
cat_pre_wt_log_scale         0.0839  1     0.7721
var_trt:cat_pre_wt_log_scale 0.2307  1     0.6310
beta:cat_pre_wt_log_scale    0.3659  2     0.8328
eclosure_m <- glmmTMB(
  eclosed ~ 
    (var_trt + beta) + cat_pre_wt_log_scale + 
    (1|session_id),
  data = d %>% 
    filter(var_trt != "constant"), 
  family = binomial(link = "logit"),
); summary(eclosure_m)
 Family: binomial  ( logit )
Formula:          
eclosed ~ (var_trt + beta) + cat_pre_wt_log_scale + (1 | session_id)
Data: d %>% filter(var_trt != "constant")

     AIC      BIC   logLik deviance df.resid 
   179.9    197.0    -84.0    167.9      122 

Random effects:

Conditional model:
 Groups     Name        Variance Std.Dev.
 session_id (Intercept) 0.1676   0.4094  
Number of obs: 128, groups:  session_id, 5

Conditional model:
                      Estimate Std. Error z value Pr(>|z|)
(Intercept)           0.551286   0.413578   1.333    0.183
var_trtlow_var        0.005064   0.373402   0.014    0.989
beta0                -0.237692   0.460689  -0.516    0.606
beta5                -0.523521   0.456587  -1.147    0.252
cat_pre_wt_log_scale  0.330169   0.220049   1.500    0.134
car::Anova(eclosure_m, type = "II")
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: eclosed
                      Chisq Df Pr(>Chisq)
var_trt              0.0002  1     0.9892
beta                 1.3212  2     0.5165
cat_pre_wt_log_scale 2.2513  1     0.1335

2.3.2 Post-hoc contrasts

2.3.3 Misc stats

r2(eclosure_m, tolerance = 10^-10)
# R2 for Mixed Models

  Conditional R2: 0.092
     Marginal R2: 0.046

2.3.4 Check model

check_mod(eclosure_m)

2.3.5 Plot

2.4 Survival

2.4.1 Model summary

surv_m_full <- coxph(
  Surv(surv_time, observed_dead) ~ 
    (var_trt + beta) * cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) +
    strata(session_id),
  data = d %>% 
    filter(var_trt != "constant") %>% 
    mutate(
      beta = as.factor(as.character(beta))
    ),
); summary(surv_m_full)
Call:
coxph(formula = Surv(surv_time, observed_dead) ~ (var_trt + beta) * 
    cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) + strata(session_id), 
    data = d %>% filter(var_trt != "constant") %>% mutate(beta = as.factor(as.character(beta))))

  n= 128, number of events= 105 

                                       coef exp(coef) se(coef)      z Pr(>|z|)
var_trtlow_var                      -0.1115    0.8945   0.2124 -0.525 0.599727
beta0                               -0.0896    0.9143   0.2579 -0.347 0.728306
beta5                               -0.2521    0.7772   0.2605 -0.968 0.333144
cat_pre_wt_log_scale                 0.9641    2.6224   0.2556  3.772 0.000162
I(cat_pre_wt_log_scale^2)            0.2450    1.2776   0.1447  1.693 0.090408
var_trtlow_var:cat_pre_wt_log_scale -0.1626    0.8499   0.1998 -0.814 0.415661
beta0:cat_pre_wt_log_scale          -0.2230    0.8001   0.2309 -0.966 0.334053
beta5:cat_pre_wt_log_scale          -0.4789    0.6195   0.2584 -1.854 0.063807
                                       
var_trtlow_var                         
beta0                                  
beta5                                  
cat_pre_wt_log_scale                ***
I(cat_pre_wt_log_scale^2)           .  
var_trtlow_var:cat_pre_wt_log_scale    
beta0:cat_pre_wt_log_scale             
beta5:cat_pre_wt_log_scale          .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                                    exp(coef) exp(-coef) lower .95 upper .95
var_trtlow_var                         0.8945     1.1179    0.5899     1.356
beta0                                  0.9143     1.0937    0.5515     1.516
beta5                                  0.7772     1.2867    0.4665     1.295
cat_pre_wt_log_scale                   2.6224     0.3813    1.5892     4.328
I(cat_pre_wt_log_scale^2)              1.2776     0.7827    0.9621     1.696
var_trtlow_var:cat_pre_wt_log_scale    0.8499     1.1766    0.5746     1.257
beta0:cat_pre_wt_log_scale             0.8001     1.2499    0.5089     1.258
beta5:cat_pre_wt_log_scale             0.6195     1.6142    0.3734     1.028

Concordance= 0.637  (se = 0.047 )
Likelihood ratio test= 22  on 8 df,   p=0.005
Wald test            = 22.31  on 8 df,   p=0.004
Score (logrank) test = 24.57  on 8 df,   p=0.002
drop1(surv_m_full, test = "Chisq")
Single term deletions

Model:
Surv(surv_time, observed_dead) ~ (var_trt + beta) * cat_pre_wt_log_scale + 
    I(cat_pre_wt_log_scale^2) + strata(session_id)
                             Df    AIC    LRT Pr(>Chi)  
<none>                          459.50                  
I(cat_pre_wt_log_scale^2)     1 460.40   2.90  0.08841 .
strata(session_id)            0 787.87 328.37           
var_trt:cat_pre_wt_log_scale  1 458.16   0.67  0.41395  
beta:cat_pre_wt_log_scale     2 458.98   3.48  0.17532  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
surv_m <- coxph(
  Surv(surv_time, observed_dead) ~ 
    (var_trt + beta) + cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) + 
    strata(session_id),
  data = d %>% 
    filter(var_trt != "constant") %>% 
    mutate(
      beta = as.factor(as.character(beta))
    ),
); summary(surv_m)
Call:
coxph(formula = Surv(surv_time, observed_dead) ~ (var_trt + beta) + 
    cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) + strata(session_id), 
    data = d %>% filter(var_trt != "constant") %>% mutate(beta = as.factor(as.character(beta))))

  n= 128, number of events= 105 

                              coef exp(coef) se(coef)      z Pr(>|z|)    
var_trtlow_var            -0.07077   0.93167  0.20845 -0.340 0.734212    
beta0                     -0.07410   0.92858  0.25808 -0.287 0.774026    
beta5                     -0.17722   0.83759  0.25804 -0.687 0.492206    
cat_pre_wt_log_scale       0.65947   1.93376  0.18164  3.631 0.000283 ***
I(cat_pre_wt_log_scale^2)  0.29584   1.34425  0.14124  2.095 0.036205 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                          exp(coef) exp(-coef) lower .95 upper .95
var_trtlow_var               0.9317     1.0733    0.6192     1.402
beta0                        0.9286     1.0769    0.5599     1.540
beta5                        0.8376     1.1939    0.5051     1.389
cat_pre_wt_log_scale         1.9338     0.5171    1.3545     2.761
I(cat_pre_wt_log_scale^2)    1.3443     0.7439    1.0192     1.773

Concordance= 0.602  (se = 0.047 )
Likelihood ratio test= 18.08  on 5 df,   p=0.003
Wald test            = 18.29  on 5 df,   p=0.003
Score (logrank) test = 19.42  on 5 df,   p=0.002
drop1(surv_m, test = "Chisq")
Single term deletions

Model:
Surv(surv_time, observed_dead) ~ (var_trt + beta) + cat_pre_wt_log_scale + 
    I(cat_pre_wt_log_scale^2) + strata(session_id)
                          Df    AIC    LRT  Pr(>Chi)    
<none>                       457.42                     
var_trt                    1 455.53   0.12 0.7341150    
beta                       2 453.89   0.48 0.7883535    
cat_pre_wt_log_scale       1 469.05  13.63 0.0002226 ***
I(cat_pre_wt_log_scale^2)  1 459.84   4.42 0.0355376 *  
strata(session_id)         0 784.70 327.28              
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.4.2 Post-hoc contrasts

2.4.3 Misc stats

r2(surv_m, tolerance = 10^-10)
  Nagelkerke's R2: 0.135

2.4.4 Check model

cox.zph(surv_m)
                          chisq df    p
var_trt                    1.60  1 0.21
beta                       1.43  2 0.49
cat_pre_wt_log_scale       2.61  1 0.11
I(cat_pre_wt_log_scale^2)  1.05  1 0.31
GLOBAL                     5.53  5 0.35

2.4.5 Plot

3 Structural Equation Model

3.1 Sub-models

3.1.1 RGR

subm_rgr <- glmmTMB(
  RGR_scale ~ 
    cat_pre_wt_log_scale + 
    cat_pre_wt_log_scale_sq + 
    sl_mean_obs_log_scale +
    prop_explore_logit_scale *  
    var_toxic_12_scale + 
    mean_toxic_conc_scale + 
    area_herb_log_scale + 
    beta_numeric_scale + 
    (1|session_id),
  data = d2
); summary(subm_rgr)
 Family: gaussian  ( identity )
Formula:          RGR_scale ~ cat_pre_wt_log_scale + cat_pre_wt_log_scale_sq +  
    sl_mean_obs_log_scale + prop_explore_logit_scale * var_toxic_12_scale +  
    mean_toxic_conc_scale + area_herb_log_scale + beta_numeric_scale +  
    (1 | session_id)
Data: d2

     AIC      BIC   logLik deviance df.resid 
   119.5    149.0    -47.7     95.5       75 

Random effects:

Conditional model:
 Groups     Name        Variance Std.Dev.
 session_id (Intercept) 0.01245  0.1116  
 Residual               0.16744  0.4092  
Number of obs: 87, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.167 

Conditional model:
                                            Estimate Std. Error z value
(Intercept)                                 -0.02766    0.07257  -0.381
cat_pre_wt_log_scale                        -3.64523    0.41435  -8.797
cat_pre_wt_log_scale_sq                     -3.01946    0.45628  -6.618
sl_mean_obs_log_scale                       -0.13590    0.05626  -2.416
prop_explore_logit_scale                     0.02370    0.04889   0.485
var_toxic_12_scale                          -0.01949    0.05115  -0.381
mean_toxic_conc_scale                       -0.25503    0.07147  -3.568
area_herb_log_scale                          0.76522    0.09726   7.868
beta_numeric_scale                          -0.11265    0.04837  -2.329
prop_explore_logit_scale:var_toxic_12_scale  0.18687    0.04488   4.164
                                            Pr(>|z|)    
(Intercept)                                 0.703071    
cat_pre_wt_log_scale                         < 2e-16 ***
cat_pre_wt_log_scale_sq                     3.65e-11 ***
sl_mean_obs_log_scale                       0.015710 *  
prop_explore_logit_scale                    0.627788    
var_toxic_12_scale                          0.703095    
mean_toxic_conc_scale                       0.000359 ***
area_herb_log_scale                         3.61e-15 ***
beta_numeric_scale                          0.019873 *  
prop_explore_logit_scale:var_toxic_12_scale 3.13e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.1.2 Variance in toxin (12 hour)

subm_var_toxic <- glmmTMB(
  var_toxic_12_scale ~ 
    sl_mean_obs_log_scale +
    beta_numeric_scale + 
    var_high + 
    (1|session_id),
  data = d2
); summary(subm_var_toxic)
 Family: gaussian  ( identity )
Formula:          
var_toxic_12_scale ~ sl_mean_obs_log_scale + beta_numeric_scale +  
    var_high + (1 | session_id)
Data: d2

     AIC      BIC   logLik deviance df.resid 
   178.2    193.0    -83.1    166.2       81 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 6.229e-11 7.892e-06
 Residual               3.956e-01 6.290e-01
Number of obs: 87, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.396 

Conditional model:
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)           -0.75075    0.09477  -7.922 2.34e-15 ***
sl_mean_obs_log_scale  0.21267    0.07751   2.744  0.00607 ** 
beta_numeric_scale    -0.15689    0.06811  -2.303  0.02126 *  
var_high               1.56021    0.13868  11.250  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.1.3 Mean step length

subm_sl <- glmmTMB(
  sl_mean_obs_log_scale ~ 
    prop_explore_logit_scale +  
    cat_pre_wt_log_scale + 
    cat_pre_wt_log_scale_sq + 
    on_toxic_logit_scale + 
    (1|session_id),
  data = d2,
); summary(subm_sl)
 Family: gaussian  ( identity )
Formula:          
sl_mean_obs_log_scale ~ prop_explore_logit_scale + cat_pre_wt_log_scale +  
    cat_pre_wt_log_scale_sq + on_toxic_logit_scale + (1 | session_id)
Data: d2

     AIC      BIC   logLik deviance df.resid 
   223.4    240.7   -104.7    209.4       80 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 4.610e-10 2.147e-05
 Residual               6.502e-01 8.064e-01
Number of obs: 87, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.65 

Conditional model:
                         Estimate Std. Error z value Pr(>|z|)   
(Intercept)               0.06163    0.09646   0.639  0.52287   
prop_explore_logit_scale -0.26175    0.08987  -2.912  0.00359 **
cat_pre_wt_log_scale      1.37948    0.70121   1.967  0.04915 * 
cat_pre_wt_log_scale_sq   1.61095    0.75569   2.132  0.03303 * 
on_toxic_logit_scale      0.21435    0.09489   2.259  0.02389 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.1.4 Mean toxin ingested

subm_toxin_ingested <- glmmTMB(
  mean_toxic_conc_scale ~ 
    beta_numeric_scale + 
    on_toxic_logit_scale + 
    var_high + 
    (1|session_id),
  data = d2 
); summary(subm_toxin_ingested)
 Family: gaussian  ( identity )
Formula:          
mean_toxic_conc_scale ~ beta_numeric_scale + on_toxic_logit_scale +  
    var_high + (1 | session_id)
Data: d2

     AIC      BIC   logLik deviance df.resid 
   145.7    160.5    -66.8    133.7       81 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 6.461e-10 2.542e-05
 Residual               2.722e-01 5.217e-01
Number of obs: 87, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.272 

Conditional model:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)           0.18349    0.07861   2.334   0.0196 *  
beta_numeric_scale   -0.11683    0.05808  -2.012   0.0443 *  
on_toxic_logit_scale  0.28426    0.06393   4.447 8.72e-06 ***
var_high             -0.80646    0.11613  -6.945 3.79e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.1.5 Mean time on toxin

subm_on_toxic <- glmmTMB(
  on_toxic_logit_scale ~ 
    var_high + beta_numeric_scale * cat_pre_wt_log_scale + 
    (1|session_id),
  family = gaussian(),
  data = d2
); summary(subm_on_toxic)
 Family: gaussian  ( identity )
Formula:          
on_toxic_logit_scale ~ var_high + beta_numeric_scale * cat_pre_wt_log_scale +  
    (1 | session_id)
Data: d2

     AIC      BIC   logLik deviance df.resid 
   230.0    247.3   -108.0    216.0       80 

Random effects:

Conditional model:
 Groups     Name        Variance Std.Dev.
 session_id (Intercept) 0.01667  0.1291  
 Residual               0.68757  0.8292  
Number of obs: 87, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.688 

Conditional model:
                                        Estimate Std. Error z value Pr(>|z|)  
(Intercept)                              0.22509    0.14140   1.592   0.1114  
var_high                                -0.41809    0.18032  -2.319   0.0204 *
beta_numeric_scale                      -0.12409    0.09766  -1.271   0.2039  
cat_pre_wt_log_scale                    -0.19306    0.13787  -1.400   0.1614  
beta_numeric_scale:cat_pre_wt_log_scale -0.23850    0.11732  -2.033   0.0421 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.1.6 Area consumption

subm_herb <- glmmTMB(
  area_herb_log_scale ~ 
    var_high * cat_pre_wt_log_scale +
    (1|session_id),
  data = d2,
); summary(subm_herb)
 Family: gaussian  ( identity )
Formula:          area_herb_log_scale ~ var_high * cat_pre_wt_log_scale + (1 |  
    session_id)
Data: d2

     AIC      BIC   logLik deviance df.resid 
   117.9    132.7    -53.0    105.9       81 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 3.501e-11 5.917e-06
 Residual               1.978e-01 4.448e-01
Number of obs: 87, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.198 

Conditional model:
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                    0.04339    0.06963   0.623  0.53313    
var_high                       0.27072    0.10127   2.673  0.00751 ** 
cat_pre_wt_log_scale           0.55603    0.08378   6.637  3.2e-11 ***
var_high:cat_pre_wt_log_scale -0.26289    0.11558  -2.275  0.02293 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.1.7 Exploration

subm_exp <- glmmTMB(
  prop_explore_logit_scale ~ 
    var_high + beta_numeric_scale * cat_pre_wt_log_scale + 
    (1|session_id),
  family = gaussian(),
  data = d2
); summary(subm_exp)
 Family: gaussian  ( identity )
Formula:          
prop_explore_logit_scale ~ var_high + beta_numeric_scale * cat_pre_wt_log_scale +  
    (1 | session_id)
Data: d2

     AIC      BIC   logLik deviance df.resid 
   250.1    267.3   -118.0    236.1       80 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 session_id (Intercept) 3.781e-10 1.944e-05
 Residual               8.830e-01 9.397e-01
Number of obs: 87, groups:  session_id, 5

Dispersion estimate for gaussian family (sigma^2): 0.883 

Conditional model:
                                        Estimate Std. Error z value Pr(>|z|)  
(Intercept)                             -0.17231    0.14397  -1.197   0.2314  
var_high                                 0.15661    0.20355   0.769   0.4416  
beta_numeric_scale                      -0.14059    0.11040  -1.273   0.2029  
cat_pre_wt_log_scale                    -0.04358    0.12362  -0.353   0.7244  
beta_numeric_scale:cat_pre_wt_log_scale  0.25977    0.13242   1.962   0.0498 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.2 d-separation test

sem_fit <- psem(
    subm_rgr,
    subm_herb,
    subm_var_toxic,
    subm_on_toxic,
    subm_toxin_ingested,
    subm_exp,
    subm_sl, 
    var_toxic_12_scale %~~% on_toxic_logit_scale, # correlated errors
    data = d2
)
sem_summary <- suppressWarnings(suppressMessages(
  summary_psem(sem_fit, 
               no_standardize_x =  c(), 
               .progressBar = FALSE, 
               direction = "var_toxic_12_scale -> area_herb_log_scale")
))
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
sem_summary$Cstat
  Fisher.C df P.Value
1   42.469 42   0.451
knitr::kable(sem_summary$dTable)
Independ.Claim Test.Type DF Crit.Value P.Value
mean_toxic_conc_scale ~ cat_pre_wt_log_scale + … coef 87 -1.5712 0.1161
var_toxic_12_scale ~ cat_pre_wt_log_scale + … coef 87 -0.1892 0.8499
prop_explore_logit_scale ~ cat_pre_wt_log_scale_sq + … coef 87 0.8624 0.3884
mean_toxic_conc_scale ~ cat_pre_wt_log_scale_sq + … coef 87 1.8002 0.0718
area_herb_log_scale ~ cat_pre_wt_log_scale_sq + … coef 87 -0.7057 0.4804
on_toxic_logit_scale ~ cat_pre_wt_log_scale_sq + … coef 87 0.7944 0.4269
var_toxic_12_scale ~ cat_pre_wt_log_scale_sq + … coef 87 0.3081 0.7580
sl_mean_obs_log_scale ~ beta_numeric_scale + … coef 87 1.1020 0.2704
area_herb_log_scale ~ beta_numeric_scale + … coef 87 0.0200 0.9840
RGR_scale ~ var_high + … coef 87 0.6991 0.4845
sl_mean_obs_log_scale ~ var_high + … coef 87 -1.5883 0.1122
on_toxic_logit_scale ~ RGR_scale + … coef 87 -1.3664 0.1718
mean_toxic_conc_scale ~ sl_mean_obs_log_scale + … coef 87 0.7459 0.4557
area_herb_log_scale ~ sl_mean_obs_log_scale + … coef 87 0.8564 0.3918
mean_toxic_conc_scale ~ prop_explore_logit_scale + … coef 87 1.0591 0.2895
area_herb_log_scale ~ prop_explore_logit_scale + … coef 87 0.3920 0.6951
on_toxic_logit_scale ~ prop_explore_logit_scale + … coef 87 0.2361 0.8134
var_toxic_12_scale ~ prop_explore_logit_scale + … coef 87 1.1110 0.2666
area_herb_log_scale ~ mean_toxic_conc_scale + … coef 87 -0.2457 0.8059
var_toxic_12_scale ~ mean_toxic_conc_scale + … coef 87 -1.3086 0.1907
on_toxic_logit_scale ~ area_herb_log_scale + … coef 87 -0.7213 0.4707

3.3 Model summary

3.3.1 Model stats

sem_summary$ChiSq
NULL
sem_summary$AIC
       AIC     AICc  K  n
1 1264.813 1276.432 51 87
sem_summary$R2
                  Response   family     link method Marginal Conditional
1                RGR_scale gaussian identity   none     0.75        0.77
2      area_herb_log_scale gaussian identity   none     0.42          NA
3       var_toxic_12_scale gaussian identity   none     0.60          NA
4     on_toxic_logit_scale gaussian identity   none     0.19        0.21
5    mean_toxic_conc_scale gaussian identity   none     0.55          NA
6 prop_explore_logit_scale gaussian identity   none     0.06          NA
7    sl_mean_obs_log_scale gaussian identity   none     0.19          NA

3.3.2 Standardized coefficients

knitr::kable(
  cbind(sem_summary$coefficients,"arrow_size"=(abs(as.numeric(sem_summary$coefficients$Std.Estimate))*20))
)
Response Predictor Estimate Std.Error DF Crit.Value P.Value Std.Estimate star arrow_size
RGR_scale cat_pre_wt_log_scale -3.6452 0.4144 87 -8.7974 0.0000 -3.5981 *** 71.962
RGR_scale cat_pre_wt_log_scale_sq -3.0195 0.4563 87 -6.6176 0.0000 -2.7708 *** 55.416
RGR_scale sl_mean_obs_log_scale -0.1359 0.0563 87 -2.4156 0.0157 -0.1448 * 2.896
RGR_scale prop_explore_logit_scale 0.0237 0.0489 87 0.4848 0.6278 0.0274 0.548
RGR_scale var_toxic_12_scale -0.0195 0.0511 87 -0.3811 0.7031 -0.0230 0.460
RGR_scale mean_toxic_conc_scale -0.2550 0.0715 87 -3.5683 0.0004 -0.2353 *** 4.706
RGR_scale area_herb_log_scale 0.7652 0.0973 87 7.8679 0.0000 0.5338 *** 10.676
RGR_scale beta_numeric_scale -0.1126 0.0484 87 -2.3287 0.0199 -0.1337 * 2.674
RGR_scale prop_explore_logit_scale:var_toxic_12_scale 0.1869 0.0449 87 4.1640 0.0000 0.2280 *** 4.560
area_herb_log_scale var_high 0.2707 0.1013 87 2.6733 0.0075 0.2315 ** 4.630
area_herb_log_scale cat_pre_wt_log_scale 0.5560 0.0838 87 6.6371 0.0000 0.7868 *** 15.736
area_herb_log_scale var_high:cat_pre_wt_log_scale -0.2629 0.1156 87 -2.2745 0.0229 -0.2795 * 5.590
var_toxic_12_scale sl_mean_obs_log_scale 0.2127 0.0775 87 2.7437 0.0061 0.1918 ** 3.836
var_toxic_12_scale beta_numeric_scale -0.1569 0.0681 87 -2.3033 0.0213 -0.1575 * 3.150
var_toxic_12_scale var_high 1.5602 0.1387 87 11.2503 0.0000 0.7875 *** 15.750
on_toxic_logit_scale var_high -0.4181 0.1803 87 -2.3186 0.0204 -0.2235 * 4.470
on_toxic_logit_scale beta_numeric_scale -0.1241 0.0977 87 -1.2706 0.2039 -0.1320 2.640
on_toxic_logit_scale cat_pre_wt_log_scale -0.1931 0.1379 87 -1.4003 0.1614 -0.1708 3.416
on_toxic_logit_scale beta_numeric_scale:cat_pre_wt_log_scale -0.2385 0.1173 87 -2.0329 0.0421 -0.2116 * 4.232
mean_toxic_conc_scale beta_numeric_scale -0.1168 0.0581 87 -2.0116 0.0443 -0.1503 * 3.006
mean_toxic_conc_scale on_toxic_logit_scale 0.2843 0.0639 87 4.4467 0.0000 0.3438 *** 6.876
mean_toxic_conc_scale var_high -0.8065 0.1161 87 -6.9447 0.0000 -0.5214 *** 10.428
prop_explore_logit_scale var_high 0.1566 0.2035 87 0.7694 0.4416 0.0808 1.616
prop_explore_logit_scale beta_numeric_scale -0.1406 0.1104 87 -1.2734 0.2029 -0.1443 2.886
prop_explore_logit_scale cat_pre_wt_log_scale -0.0436 0.1236 87 -0.3526 0.7244 -0.0372 0.744
prop_explore_logit_scale beta_numeric_scale:cat_pre_wt_log_scale 0.2598 0.1324 87 1.9617 0.0498 0.2224 * 4.448
sl_mean_obs_log_scale prop_explore_logit_scale -0.2618 0.0899 87 -2.9124 0.0036 -0.2840 ** 5.680
sl_mean_obs_log_scale cat_pre_wt_log_scale 1.3795 0.7012 87 1.9673 0.0491 1.2778 * 25.556
sl_mean_obs_log_scale cat_pre_wt_log_scale_sq 1.6109 0.7557 87 2.1318 0.0330 1.3873 * 27.746
sl_mean_obs_log_scale on_toxic_logit_scale 0.2144 0.0949 87 2.2589 0.0239 0.2244 * 4.488
~~var_toxic_12_scale ~~on_toxic_logit_scale 0.5323 - 87 5.7626 0.0000 0.5323 *** 10.646

3.3.3 Post hoc contrast

contrast_by_pre_wt(subm_herb, "var_high", "trt")
$emmeans
cat_pre_wt_log_scale = -2:
 var_high emmean    SE df lower.CL upper.CL
        0 -1.069 0.200 81   -1.467   -0.670
        1 -0.272 0.198 81   -0.666    0.122

cat_pre_wt_log_scale =  2:
 var_high emmean    SE df lower.CL upper.CL
        0  1.155 0.161 81    0.836    1.475
        1  0.900 0.150 81    0.603    1.198

Confidence level used: 0.95 

$contrasts
cat_pre_wt_log_scale = -2:
 contrast              estimate    SE df t.ratio p.value
 var_high1 - var_high0    0.797 0.281 81   2.830  0.0059

cat_pre_wt_log_scale =  2:
 contrast              estimate    SE df t.ratio p.value
 var_high1 - var_high0   -0.255 0.219 81  -1.162  0.2485
sd(subm_herb$frame$var_high)/sd(subm_herb$frame$area_herb_log_scale)
[1] 0.8549509
emmeans::emtrends(subm_exp, 
                 ~ beta_numeric_scale | cat_pre_wt_log_scale, 
                 var = "beta_numeric_scale", 
                 at = list(cat_pre_wt_log_scale = c(-2,2)))
cat_pre_wt_log_scale = -2:
 beta_numeric_scale beta_numeric_scale.trend    SE df lower.CL upper.CL
           7.88e-17                   -0.660 0.323 80    -1.30  -0.0173

cat_pre_wt_log_scale =  2:
 beta_numeric_scale beta_numeric_scale.trend    SE df lower.CL upper.CL
           7.88e-17                    0.379 0.246 80    -0.11   0.8677

Results are averaged over the levels of: var_high 
Confidence level used: 0.95 
sd(subm_exp$frame$beta_numeric_scale) * sd(subm_exp$frame$prop_explore_logit_scale)
[1] 0.974435
emmeans::emtrends(subm_on_toxic, 
                 ~ beta_numeric_scale | cat_pre_wt_log_scale, 
                 var = "beta_numeric_scale", 
                 at = list(cat_pre_wt_log_scale = c(-2,2)))
cat_pre_wt_log_scale = -2:
 beta_numeric_scale beta_numeric_scale.trend    SE df lower.CL upper.CL
           7.88e-17                    0.353 0.286 80   -0.217    0.923

cat_pre_wt_log_scale =  2:
 beta_numeric_scale beta_numeric_scale.trend    SE df lower.CL upper.CL
           7.88e-17                   -0.601 0.217 80   -1.034   -0.169

Results are averaged over the levels of: var_high 
Confidence level used: 0.95 
sd(subm_on_toxic$frame$beta_numeric_scale) * sd(subm_on_toxic$frame$on_toxic_logit_scale)
[1] 0.9403361

3.4 SEM plots

knitr::include_graphics(paste0(getwd(), "/graphs/manuscript1_figures/SEM_simplified.png"), rel_path = FALSE)

4 Session Info

sessionInfo()
R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)

Matrix products: default


locale:
[1] LC_COLLATE=English_United States.utf8 
[2] LC_CTYPE=English_United States.utf8   
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

time zone: America/New_York
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] spat1f_0.0.1       herbivar_0.2.1     forcats_0.5.1      stringr_1.5.0     
 [5] dplyr_1.1.2        purrr_1.0.1        readr_2.1.2        tidyr_1.3.0       
 [9] tibble_3.2.1       tidyverse_1.3.1    extraDistr_1.9.1   imager_0.42.13    
[13] magrittr_2.0.3     doSNOW_1.0.20      snow_0.4-4         iterators_1.0.14  
[17] foreach_1.5.2      glmmTMB_1.1.7      tidybayes_3.0.2    performance_0.10.4
[21] sjPlot_2.8.15      ggpubr_0.6.0       ggplot2_3.4.2      piecewiseSEM_2.3.0
[25] survival_3.3-1    

loaded via a namespace (and not attached):
  [1] svUnit_1.0.6             splines_4.3.1            later_1.3.0             
  [4] bitops_1.0-7             cellranger_1.1.0         polyclip_1.10-0         
  [7] reprex_2.0.1             lifecycle_1.0.3          sf_1.0-8                
 [10] Rdpack_2.3.1             rstatix_0.7.2            doParallel_1.0.17       
 [13] rprojroot_2.0.3          vroom_1.5.7              processx_3.6.1          
 [16] lattice_0.21-8           MASS_7.3-60              insight_0.19.3          
 [19] readbitmap_0.1.5         ggdist_3.1.1             backports_1.4.1         
 [22] sass_0.4.1               rmarkdown_2.24           jquerylib_0.1.4         
 [25] yaml_2.3.5               remotes_2.4.2            httpuv_1.6.5            
 [28] gap_1.2.3-6              sp_1.5-0                 sessioninfo_1.2.2       
 [31] pkgbuild_1.3.1           cowplot_1.1.1            DBI_1.1.3               
 [34] minqa_1.2.4              RColorBrewer_1.1-3       lubridate_1.8.0         
 [37] DHARMa_0.4.6             multcomp_1.4-25          abind_1.4-5             
 [40] pkgload_1.3.2            rvest_1.0.2              imagerExtra_2.0.0       
 [43] TH.data_1.1-1            tensorA_0.36.2           sandwich_3.0-2          
 [46] spatstat.utils_3.0-1     units_0.8-0              fitdistrplus_1.1-11     
 [49] codetools_0.2-19         MuMIn_1.47.5             xml2_1.3.3              
 [52] tidyselect_1.2.0         ggeffects_1.2.3          farver_2.1.0            
 [55] lme4_1.1-29              matrixStats_1.2.0        stats4_4.3.1            
 [58] jsonlite_1.8.8           e1071_1.7-11             ellipsis_0.3.2          
 [61] bmp_0.3                  emmeans_1.7.5            tools_4.3.1             
 [64] Rcpp_1.0.8.3             glue_1.6.2               mgcv_1.8-40             
 [67] xfun_0.39                distributional_0.3.0     usethis_2.1.6           
 [70] withr_2.5.0              numDeriv_2016.8-1.1      fastmap_1.1.0           
 [73] boot_1.3-28.1            fansi_1.0.3              callr_3.7.0             
 [76] digest_0.6.29            R6_2.5.1                 mime_0.12               
 [79] estimability_1.3         colorspace_2.0-3         spatstat.data_3.0-0     
 [82] jpeg_0.1-9               see_0.7.1                DiagrammeR_1.0.9        
 [85] utf8_1.2.2               generics_0.1.2           class_7.3-22            
 [88] prettyunits_1.1.1        httr_1.4.3               htmlwidgets_1.5.4       
 [91] pkgconfig_2.0.3          gtable_0.3.0             circular_0.5-0          
 [94] moveHMM_1.9              htmltools_0.5.2          carData_3.0-5           
 [97] profvis_0.3.7            fftwtools_0.9-11         TMB_1.9.4               
[100] scales_1.2.1             png_0.1-7                posterior_1.2.2         
[103] CircStats_0.2-6          knitr_1.43               rstudioapi_0.13         
[106] geosphere_1.5-18         tzdb_0.3.0               coda_0.19-4             
[109] visNetwork_2.1.0         checkmate_2.1.0          nlme_3.1-162            
[112] nloptr_2.0.3             proxy_0.4-27             cachem_1.0.6            
[115] zoo_1.8-12               sjlabelled_1.2.0         KernSmooth_2.23-20      
[118] parallel_4.3.1           miniUI_0.1.1.1           desc_1.4.1              
[121] pillar_1.9.0             grid_4.3.1               vctrs_0.6.3             
[124] urlchecker_1.0.1         promises_1.2.0.1         car_3.1-2               
[127] dbplyr_2.2.0             arrayhelpers_1.1-0       xtable_1.8-4            
[130] evaluate_0.15            amt_0.2.1.0              mvtnorm_1.1-3           
[133] cli_3.6.1                compiler_4.3.1           rlang_1.1.1             
[136] crayon_1.5.1             ggsignif_0.6.3           labeling_0.4.2          
[139] ggmap_3.0.2              modelr_0.1.8             classInt_0.4-7          
[142] RcppArmadillo_0.11.2.0.0 ps_1.7.1                 plyr_1.8.7              
[145] sjmisc_2.8.9             fs_1.5.2                 gap.datasets_0.0.5      
[148] stringi_1.7.12           deldir_1.0-6             assertthat_0.2.1        
[151] munsell_0.5.0            tiff_0.1-11              spatstat.geom_3.0-3     
[154] devtools_2.4.5           bayestestR_0.13.1        Matrix_1.5-4.1          
[157] sjstats_0.18.2           hms_1.1.1                bit64_4.0.5             
[160] qgam_1.3.4               RgoogleMaps_1.4.5.3      shiny_1.7.1             
[163] highr_0.9                haven_2.5.0              rbibutils_2.2.8         
[166] igraph_1.3.2             broom_1.0.5              memoise_2.0.1           
[169] bslib_0.3.1              bit_4.0.4                readxl_1.4.0